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We study, in random sparse networks, finite size scaling of the spin glass susceptibility xsg, 
which is a proper measure of the de Almeida-Thouless (AT) instability of spin glass systems. Using a 
phenomenological argument regarding the band edge behavior of the Hessian eigenvalue distribution, 
we discuss how xsg is evaluated in infinitely large random sparse networks, which are usually 
identified with Bethe trees, and how it should be corrected in finite systems. In the high temperature 
region, data of extensive numerical experiments are generally in good agreement with the theoretical 
values of xsg determined from the Bethe tree. In the absence of external fields, the data also show a 
scaling relation xsg = N 1 ^ 3 F(N 1,la \T — T c \/T c ), which has been conjectured in the literature, where 
T c is the critical temperature. In the presence of external fields, on the other hand, the numerical 
data are not consistent with this scaling relation. A numerical analysis of Hessian eigenvalues implies 
that strong finite size corrections of the lower band edge of the eigenvalue distribution, which seem 
relevant only in the presence of the fields, are a major source of inconsistency. This may be related 
to the known difficulty in using only numerical methods to detect the AT instability. 

PACS numbers: 75.50.Lk, 75.40Mg, 64.60.De 



I. INTRODUCTION 



The discovery of spontaneous replica symmetry break- 
ing (RSB) in the low temperature region of spin glass 
(SG) models^ is perhaps the most significant achievement 
in the field of statistical mechanics of disordered systems. 
It was first found in the analysis of the Sherrington- 
Kirkpatrick (SK) model, which is a fully connected mean 
field SG model, by introducing the Parisi ansatz for de- 
scribing an RSB state as a relevant saddle point of the 
replicated free energy, which is intrinsically replica sym- 
metric. The mathematical technicalities of construct- 
ing the solution generated tremendous controversy in 
the early days of SG research. However, alternative 
mean field approaches^ and mathematically rigorous 



arguments^ now support the correctness of the Parisi 
solution. These days, it is commonly accepted that RSB 
does occur in a class of mean field models although the 
existence of RSB in systems of finite dimension is still 
under debate^. 

In the present situation, SG models on random sparse 
networks play a special role. A random sparse network 
is constructed such that each node is randomly coupled 
to a finite number of other nodes. As the random con- 
struction guarantees statistical uniformity of the network 
structure, the random sparse model is classified as a mean 
field model that exhibits RSB. At the same time, unlike 
fully connected models, a concept of adjacency between 
nodes is naturally introduced in the random sparse net- 
work, which may make it possible to characterize the 
RSB transition by using the concept of correlation length 
(as in the standard analysis of models in finite dimen- 
sions). For very large (possibly infinitely large) random 
sparse networks, a solution for the SG model exists when- 



ever the SG correlation decays fast enough (in practice, in 
the replica symmetric phases). This solution is obtained 
by assuming equivalence between the random sparse net- 
work and the corresponding Bethe tree. 

In SG models that undergo a continuous phase tran- 
sition as the temperature decreases, the onset of RSB is 
signaled by the de Almeida-Thouless (AT) instability^, 
which corresponds to the divergence of the SG suscep- 
tibility. The AT instability naturally defines a critical 
line T C (H) in the temperature (T) - external field (H) 
plane. The shape of this line in the fully connected SK 
model is very peculiar: H C (T) diverges for T — > 0; that 
is, an RSB phase exists for any value of the external field. 
This behavior is rather unusual and can take place only 
in models where the coordination number diverges in the 
thermodynamic limit. In more realistic situations, T C (H) 
should become when the external field reaches a critical 
value H c : SG models on random sparse networks display 
this more realistic behavior. 

A very interesting (and still largely open) question is 
how the AT instability is observed in finite systems. Nu- 
merical studies of the SK model have failed to identify the 
critical point with H ^ and have reported very strong 
finite size effects^. For SG models on random sparse net- 
works, much less is known: there are very few numerical 
studies^ and a validation of the finite size scaling relations 
is still lacking. 

In this paper, we mainly examine how the critical con- 
dition for the AT instability obtained in the thermody- 
namic limit should be corrected for finite sized systems. 

A proper measure for detecting the AT instability is 
the spin glass susceptibility defined as2£ 

i,3 
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where N is the number of spins, Si is a spin variable, and 
(• • • ) (respectively, 7TT ) denotes thermal (respectively, 
configurational or disorder) averages. For T > T C (H), 
in the infinitely large system limit, a random sparse net- 
work can be accurately approximated by a Bethe tree, 
thus providing a direct expression of XSG- On the other 
hand, for large but finite systems, that expression for xsg 
has to be corrected in an appropriate manner to account 
for finite size effects. We will show that the edge behav- 
ior of the eigenvalue distribution p(A) of the susceptibility 
matrix plays a key role in this correction. That is, when 
p(X) oc (A — A m i n ) 7 holds in the vicinity of the lower band 
edge of the distribution, i.e., near A m i n , the finite size cor- 
rection amounts to a scaling relation xsg = N^F^N 1 "), 
where u = (1 - + 7) and t = (T - T c )/T c , unless 

Amin is significantly influenced by finite size effects. In the 
absence of external fields, the results of random matrix 
theory indicate 7 = 1/2, which yields a known relation 
for mean field SG models xsg = N^FitN 1 / 3 ). Data 
obtained from extensive numerical experiments confirm 
that this relation is fairly accurate. On the other hand, 
numerical data for N < 2 10 show a considerable discrep- 
ancy with the same scaling relation in the presence of 
external fields even when ui is optimally tuned. Numeri- 
cal analysis on the Hessian for relatively smaller systems 
of N < 2 s indicates that A m i n has strong finite size cor- 
rections in the presence of fields, whereas the profile near 
Amin of the eigenvalue distribution does not change much, 
which makes it practically difficult to identify T c by us- 
ing only numerical methods. This may be the reason why 
the AT instability is hard to observe in finite-dimensional 
models. 

This paper is organized as follows. The next section in- 
troduces the model to be examined. In section 3, we show 
how xsg is evaluated for infinitely large random sparse 
networks. We also derive the finite size scaling relation 
on the basis of phenomenological considerations about 
the eigenvalue distribution of the susceptibility matrix. 
In section 4, we discuss the numerical experiments ex- 
amining the validity of the results obtained in section 3. 
The final section is devoted to a summary. 



II. MODEL DEFINITION 

We will study SG models defined on C-regular ran- 
dom graphs in the presence of an external field H. The 
Hamiltonian is given by 



ns) = - £ 
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where E is the set of edges in the random graph, which 
is chosen uniformly among all graphs of N nodes and 
M = NC/2 edges, having exactly C edges per node. 
Each node contains an Ising spin Si G {+1,-1} (i — 
1, . . . ,N), and the couplings arc quenched i.i.d. random 



variables extracted from 



P(J) = -[S(J 



1) + 5(J + 1) 



(2) 



As already mentioned, we denote the thermal 
averages with respect to the canonical distri- 
bution of inverse temperature /3 = T _1 as 
(•••) = £ 5 (---)exp[-/3ft(S)]/Z(/3), where 
— Es cx p[ _ ^(^)] i s the partition func- 
tion. Configurational averages with respect to the 
generation of the couplings and of the graph are denoted 
as (•••)■ 



III. DE ALMEIDA-THOULESS INSTABILITY 

A. Infinitely large systems 

Thermal averages arc difficult to evaluate computa- 
tionally. However, when a given graph is free of cycles 
(i.e., it is tree), it is known that the Bethe approxima- 
tion is exact and the belief propagation (BP) algorithm 
provides exact averages in a practical time^l. The BP 
iterative method is as follows: 



tanh 



tanh(/3 J i3 ) tanh fiH 
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where 



(3) 

represent message variables which are trans- 
mitted between nodes, di is the set of neighbors of i, 
and \j stands for exclusion of j. From the fixed point of 
Eq. ([3]), the thermal average of spin Si is evaluated as 



mi = (Si) = tanh 



/3H 



E 1 



(4) 



The results obtained by the BP algorithm, Eq. (|3|), 
are just an approximation for general lattices with loops. 
However, the length of cycles in a random sparse network 
typically grows as O(lnTV) as the size of the graph N 
tends to infinity. This implies that for high temperatures, 
T > T c , where the spatial correlations decay fast enough, 
thermal averages in sufficiently large random networks 
can be accurately evaluated by using BP, i.e. as if they 
were computed on a tree. This idea is also useful for 
analyzing the AT instability. 

Let us focus on a spin Si in a large random sparse 
network and approximate the lattice with a tree rooted 
at i. A distinctive property of the tree is that any pair 
of points is linked by a unique path. Noticing this, we 
can define a path connecting i with j, placed at distance 
G from i, and assign a label g = 0,1, . . . ,G to the nodes 
along the path from i to j (g = and g = G correspond 
to i and j, respectively). On this tree, the two-point 
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correlation between i and j can be exactly computed as 

_ x drrii 



{S i S j )-(S i )(S j ) = 
= 
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(5) 



Here, the derivative with respect to the field Hj acting 
on j has been replaced by the one with respect to any BP 
message arriving at site j. Statistical uniformity guaran- 
tees that the configurational average of the square of Eq. 
© 



({SiSj-iSiHSi))* 



drrii 
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depends only on G. On the other hand, the number of 
spins at a distance G from % is given by C(C — 1) G_1 on 
the tree. Thus, the spin glass susceptibility xsg can be 
evaluated as 



Xsg = (1 - mf) 



J2 C(C - lf-\{S S G ) - (So) (S G )f 



G-l 



l-(C-l)e-*' 
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for T approaching T c from above. The quantity ^ is the 
inverse of the SG correlation length 
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From Eq. © , we obtain the condition for the divergence 
of Xsg 



(C-l)e- 



1 



(8) 



giving the AT instability for infinitely large system s 13 ' 14 . 

Three points are noteworthy. First, unlike the usual 
critical phenomena in finite dimensions, the AT instabil- 
ity on random sparse networks is not accompanied by 
the divergence of the correlation length: even at the crit- 
ical point, the correlation length is finite and equal to 
£ = "f^ 1 = l/ln(C — f). This is because the num- 
ber of nodes at a distance G from a fixed node in a 
random network grows exponentially fast and propor- 
tionally to (C — I) G , which is not the case in models 
of finite dimensions. Second, for H = 0, we can ob- 
tain an analytical expression for xsg- Indeed, as long 
as (C — I) tanh(/3) 2 < I, a unique convergent solution of 



H->J 



0, which in turn 



Eq. (|3|) exists and is given by 
implies ((S S G ) - (S Q ) (S G )f = [tanh(^)] 2G and 



XSG 



I + tanh(/?) 5 



I - (C- l)tanh(/3) 2 
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The above equation yields the critical temperature by 
setting (C — I)tanh(/3 C ) 2 = I, and it shows that as 
T — > T c from above, xsg diverges as 0(|t| -1 ), where 
t = (T — T c )/T c . in perfect agreement with the known 
AT instability in the absence of an external fiel d 15 ' 16 . 
Third, one can still numerically assess Eq. © in a prac- 
tical time even in the presence of an external field. For 
this, we utilize a property of BP operating on a Bethc 
tree whereby the distributions of the message variables 
for typical sample systems can be obtained as a set of 
solutions of functional equations, 



7r(w) 
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(10) 

where f(x,y) = tanh ( tanh(x) tanh(y)). Eq. (fTU)) can 
be solved numerically with a sampling method in a rea- 
sonable time^. This implies that, for given G, a sample 
of (SqS g ) — (Sq) (S G ) can be generated by Eq. {5) from 
the following BP dynamics: 



l g->g-i 



f{PJg, Ug+l^-g+Vg) 



Iff) 



where g = G, G — 1, . . . , 1, and J g and r g represent in- 
dependent random numbers respectively sampled from 
Eq. ^ and from the distribution 

.C-2 / C-2 \ 

Kr)= / l[du M ir(u M )slr-PH-^2uA . (12) 

Here, r g stands for the sum of messages from C — 2 
branches that merge with the <?-th node on the path. 
The computational cost of this evaluation scales as 
0(G) per sample, which is computationally feasible. 
Therefore, for fixed G, one can numerically evaluate 

((SqS g ) — (Sq) (S G )) 2 by using the sampling stochastic 
process of Eqs. (|TT1) (|I2[) many times. Estimates for 
several G can be then interpolated with an appropriate 
polynomial in 1/G. This makes it possible to practi- 
cally compute Eq. ([7]) by extrapolation of 1/G —> in 
the fitted polynomial. Table Q] lists estimates of T c for 
H = 0, 0.1 , 0.2 and 0.3 for the case of C = 4. The values 
for H were evaluated by extrapolating fourth degree 
polynomials fitted to data of G = 1, 2, . . . , 20 that were 
computed from 10 7 samplings. T c = 1.5187 . . . for H = 
can be obtained by solving 3tanh 2 (I/T c ) = 1. \& varies 
linearly with respect to t around the critical temperature 
T c . This implies that xsg scales as 0(i -1 ) close to T c in 
the limit of N -> oo. 

We also checked that the method for computing the 
AT line is equivalent to the method based on perturb- 
ing the BP messages and then observing the subsequent 
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H 



0.0 



0.1 



0.2 



0.3 



1.51865 1.3053 1.1808 1.0770 



TABLE I: The critical temperatures T c for fixed connectivity 
4 and four different external fields H . 



evolution of the perturbation, the critical temperature 
being defined as the lowest temperature such that the 
perturbation does not grow under BP iteration^. 

Note that although the graph is regular, i.e. all the ver- 
tices are equivalent, the presence of a uniform field pro- 
duces many heterogeneities and does not allow for the ex- 
istence of a factorized solution (and this may eventually 
lead to larger fluctuations in finite systems^) . The main 
effect is that the two-spin correlation (SiSj) — (Si)(Sj) 
fluctuates a lot between different pairs of spins separated 
by the same distance. Figure Q] plots the rate function 
Q,(v) = G^ 1 \\\P(v) of the probability distribution P(v) 
for the logarithm of the correlation 



v = lim — In 

G^oo G 



((SoS G ) - (S )(S G )Y 



Correlations contributing the most to the assessment of 



= — argmax 



are marked by dots in Fig. [TJ 

and are clearly larger than the most probable correlations 
corresponding to the maximum of Q(v). This means that, 
as soon as H ^ 0, the long-range order in the model is 
produced essentially by very few pairs of strongly cor- 
related spins, while the vast majority of pairs of spins 
remain uncorrelated. 
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FIG. 1: Profiles of rate function at the AT critical- 

ity for several values of external field H in the case of 
C = 4. Values of the critical temperatures are shown 
in table [T] In order to evaluate Q(u), we first computed 

(a) = - lime-Wl/G) In {(S S G ) -(S )(S G )) 23 by extrap- 
olating numerical data for G = 1, 2, . . . , 20 to G — > oo. Ap- 
plying a Legendre transformation to this function yields 
as follows: v = —(d/ds)^(s) and fl — —sv — 9(s), where Q 
has been parameterized by a conjugate variable s. For draw- 
ing the profiles shown in the figure, we numerically evaluated 

{(SoS G ) - (So) {S G )) 2a based on 10 7 samples of eq. (HQ and 
varied s in the range of — 1 < s < 9. The profiles for H ^ 
indicate that the dominant values of v for the AT criticality 
(dots) are considerably larger than the most probable values 
of v (crosses). The physical implication of this is that the 
AT instability for H 7^ is induced by a small number of 
atypically large spin correlations. 



B. Finite systems 

So far, we have reviewed how T c for the AT instability 
can be evaluated by computing xsg m infinitely large 
random sparse networks. However, xsg is intrinsically 
upper-bounded by N and never diverges as long as N is 
finite. Therefore, we have to examine how the scenario in 
the previous section should be modified in finite systems 
so that we can appropriately analyze data from numerical 
experiments. 

A naive correction taking the finiteness of the system 
into account can be made by truncating the summation 
of Eq. © at a finite number G = G max , which leads to 



XSG oc 



l _ (( C_l )e -* ) 
1 - (C- l)e 



(13) 



since G cannot tend to infinity when the system is fi- 
nite. Unfortunately, the following considerations indi- 
cate that such a correction is not appropriate for describ- 
ing the behavior in the vicinity of T c . The right-hand 
side of Eq. (fTS]) gives G max when the critical condition 
1 — (C — l)e~* — » holds. G max grows monotonically as 
increases. However, the growth rate is only 0(lnN) 
since N ~ C(C — l) Gmax must hold for satisfying the 



constraint concerning the number of nodes. This rate 
is obviously too slow since numerical experiments show 
that the xsg of finite systems grows as 0(N 1 ^ 3 ) at the 
critical condition, at least, for H = 0^. 

This discrepancy indicates that effects of self-inter- 
actions, which are ignored in the Bethe tree approxima- 
tion, must be taken into account when evaluating the 
dependence of xsg on the system size N in the vicinity 
of T c . Unfortunately, such an evaluation requires a com- 
plicated calculation and still does not lead to an accurate 
expression in general. Therefore, to avoid these technical 
difficulties, we shall employ a phenomenological deriva- 
tion. 

Consider the Hessian A = x •> where x denotes 
a susceptibility matrix with elements Xij — i(SiSj) ~ 
(Si) (Sj)). As a working hypothesis, we assume that the 
eigenvalues of A, Ai < A2 < . . . < Ajv, obey a continuous 
distribution p(A), which behaves as 



p(X) oc (A - A min ) n 



(14) 



close to the lower band edge A m i n for T > T c and N — > 00. 
Moreover, we shall assume for the moment that A m i n is 
not heavily modified by finite corrections. For H = 0, 
an analysis of random matrices of fixed weights, in con- 
junction with Thouless-Anderson-Palmer theory^, im- 
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plies that the distribution can be expressed as 



p(A; A,/x) 



1 V4(C - 1)A 2 - (A - M )s 
2tt CA 2 - (A-/i) 2 /C 



(15) 



using certain parameters A and which supports 7 = 
1/2 for C > 3. However, here, we do not exclude the 
possibility that 7 may depend on T and H. Equation (fT4")l 
provides another expression of xsg 



XSG 



1 



N 



fe=i 



dAA- 2 p(A)ocA-£-^ (16) 



as AT — > 00. Assuming that xsg diverges as 0(|i| _1 ) at 
criticality, A m i n oc £ 1 /( 1_ t) holds as T approaches T c from 
above in the limit of A" — > 00. 

However, the statistical fluctuations of the eigenvalues 
are not negligible around T c for large but finite N. As 
a first approximation, therefore, let us regard Afc (k = 
1,2, ... , N) as independently and identically distributed 
(i.i.d.) random variables extracted from p(A). Since Ai is 
the smallest value among the A" i.i.d. random variables, 
the theory of extreme value statistics^ indicates that 
magnitude of the fluctuation of Ai can be evaluated by a 
simple equation: 



A 



Ai 



dAp(A)~0(l), 



(17) 



which yields a scaling relation Ai — A m i n oc A" 1/(1+7). 

Replacing A m j n by its scaling relation in terms of 
tV(i-7) ; Eq. (TTT]) leads to the following expression for 
the smallest eigenvalue: 

Ai = A 1 t 1 /(i-7) +A r-V(i+7)^ 1 = 

= ^-1/(1+7) ^(tJV^jVCi-T) +^ ; (is) 

where A\ is a constant and £1 is a random variable taking 
values O(l). This derivation also indicates that Afc can 
be expressed similarly to Eq. (fT5|) as long as k ~ O(l). 
Accordingly, all contributions to xsg from Afc with k ~ 
O(l) can be summed together in a scaling relation like 



1 

N 



E 

fc~0(l) 



N^g ItN^ 



after being averaged with respect to the £fc. Here, g(x) is 
a well-behaved function which returns O(l) constant for 
x = and decays polynomially as 1/x for x S> 1. 

The contribution from all larger eigenvalues, Afc with 
k O(N), to XSG can be written in an integral form 
similar to Eq. (|16|) by substituting the lower band edge 
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l\ nliD +0(N- 1 /( 1 +-')) 

oc (A min + 0(Ar-v(i+7)^ 
- r x h (tN^^j . (19) 



In the last relation, we have used again the scaling rela- 
tion for A m ; n , and h(x) is another well-behaved function 
that is proportional to x for \x\ <C 1 and converges to a 
certain constant as x —> 00. Combining the two contri- 
butions, we obtain the finite size scaling relation of %SG 
as 



Xsg - N^g (tN^J +t~ 1 h (tN^+$^ 



1 — 1 / 1 — 1 
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where F(x) = g(x) + h(x)/x. The properties of g(x) 
and h(x) guarantee that F(0) is a finite constant and 
F(x) - 1 + 0{x~ 1 ) for x > 1. 

For H = and C > 3, 7 = 1/2 yields the scal- 
ing law xsg = N 1 / 3 F(tN 1 / 3 ). This relation was often 
assumed in earlier studies on SG models of the mean 
field type^ i 21 ' 22 . However, as far as the authors know, 
there has been no numerical validation of this relation, 
in particular, for the scaling exponent with respect to 
t = \T - T c \/T c , even for the case of H = 0. In ad- 
dition, there is no theoretical guarantee that 7 = 1/2 
always holds for H > case. As there are only few an- 
alytical schemes available for dealing with SG models of 
finite dimension, we need to build up a solid basis for nu- 
merical studies. Circumstantially comparing the results 
of numerical experiments and theoretical predictions of 
Eqs. ([5]) and (|2U)) for the current system is a great step 
towards fulfilling this purpose. 



IV. NUMERICAL EXPERIMENTS 

In order to verify the above-mentioned behavior 
around the AT instability, we performed large numer- 
ical experiments on systems with C — 4 and sizes 
A^ = 2 5 , 2 6 , . . . , 2 10 , using the replica exchange (par- 
allel tempering) Markov chain Monte Carlo (MCMC) 
metho d 23 ' 24 . Apart from some test runs on small sys- 
tems with H — 0, we ran extensive simulations on fields 
H = 0.1,0.2,0.3 at respectively 33, 34, and 36 differ- 
ent temperatures distributed around T c . For equilibrat- 
ing the systems, we performed 2 21 MC sweeps (MCSs) 
and computed thermal averages from 2 21 more MCSs af- 
ter the equilibration time. Equilibration was tested by 
comparing the averages obtained by using half and one 
quarter of the total MCSs. To accelerate equilibration, 
replicas of adjacent temperatures were exchanged once 
every 30 MCSs. We simulated 16000 samples for each 
size. 

Figure [5] shows the results of runs with H = 0. The 
spin glass susceptibility rescaled by a factor N^ 1 / 3 (corre- 
sponding to 7 = 1/2) nicely crosses at the critical temper- 
ature, as predicted analytically. The inset should show 
the scaling function (if finite size effects were absent), but 
we can clearly see that the data collapse is good only in 
the high temperature (low j3 = T _1 ) region. 

Next, let us turn to the case of external fields. Ex- 
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FIG. 2: Rescaled spin glass susceptibility xsg without exter- 
nal field. 



panding ((S^) - (Si) (S,)) 2 as 

(SiSj) (SiSj) - 2 (SiSj) (Si) (Sj) + (Si) (Sj) (Si) (S^ 

and using different real replicas for computing different 
thermal averages at the same time, the above equation 
becomes 

(sl Sj Sf Sj) — 2 (Si Si Si Si) + (SlSjSfS*) , 
we can write the spin glass susceptibility xsg &s 



XSG = N((qf 2 ) - 2 (912913) + (913924)) 



N l(q 2 12 ) - 2 (q 12 q 13 ) + (g 12 )" 



(21) 



where q a b = N^ 1 J2i=i ^i^i 1S the overlap between two 
replicas. Actually we computed xsg m Eq.()21j) by mea- 
suring overlaps from four different replicas a, b = 1, 2, 3, 4 
evolving independently, in order to reduce correlations ef- 
fects and the noise-to-signal ratio. 

Figure [3] compares the susceptibilities xsg measured 
numerically on systems of sizes N = 2 6 ,2 8 ,2 10 with 
the ones computed analytically on the Bethe tree for 
H = 0.1,0.2, and 0.3. We see that the numerical data 
at high temperatures converge nicely to the theoreti- 
cal estimates on the Bethe tree. Note that the diam- 
eter of the regular random graph with C = 4 is only 
ln(JV)/ln(C - 1) ~ 6.3 even for the case of N = 2 lri . 
This indicates that accuracy of the Bethe approxima- 
tion is not determined only by the size of the graphs 
or, more precisely, by the length of the shortest loops. 
The relative strength of the self-interactions compared 
to the size of the graphs plays a key role in determining 
the accuracy. In other words, even if the graph contains 
many loops, which may significantly contribute to self- 
interaction terms (that are missing on trees) , the lack of 
correlation in the topology makes the net contribution 
of these loops very small. The final result is that the 
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FIG. 3: Inverse of the spin glass susceptibility versus temper- 
ature for fields H = 0.1 (top), H = 0.2 (middle), and H = 0.3 
(bottom). 



critical window size for the susceptibilities scales as an 
inverse power of N rather than 1/ In N. 

Figure [3] (top) shows the analytical curves correspond- 
ing to H = and H = 0.1. One can see how the 
H = 0.1 data closely follow the H = curve as long 
as T > T C {H = 0). Only below T C (H = 0) does the 
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data change its curvature and acquire the correct linear 
behavior in T — T C (H = 0.1). Unfortunately, this change 
happens at very large values of xsg- an d thus, the asymp- 
totic scaling behavior may be difficult to observe. For 
larger fields, H = 0.2 and H = 0.3, the influence of the 
H = fixed point is weaker and the theoretical suscepti- 
bility curves are qualitatively similar to the H = curve, 
with the linear part in T — T c extending over a wider 
range. Nonetheless, for these larger fields, the values of 
Xsg are much smaller, and thus, the asymptotic behavior 
may be difficult to observe in this case as well. 

We tried to identify the critical point, T c . from the 
finite size scaling of the numerical data, by using 7 = 
1/2 for any field value. For this, we plotted N~ 1 / 3 xsg 
versus temperature and looked for a crossing point of 
data sets having different TV, which should correspond 
to T c in the thermodynamic limit. We see from Fig. [4] 
that finite size corrections are rather large, especially for 
H = 0.1 and H = 0.3, and change sign depending on 
the value of the field. The insets in Fig. [4] zoom in on 
the region containing the crossings for all data, to reveal 
whether or not the crossings move towards the analytical 
T c value computed under the tree approximation. The 
H = 0.1 crossing points move in the right direction, but 
they do so very slowly; most probably due to the H = 
fixed point at T C (H = 0) ~ 1.52 in the vicinity. The 
H = 0.3 crossing points also move towards T c and do so 
faster than those of H = 0.1, although they come from 
the low temperature phase. The H = 0.2 crossing points 
are more complex, because the crossing point apparently 
move leftward, away from the critical point. The most 
natural explanation for this behavior is that for larger 
system sizes the crossing point moves back again towards 
T c , as in the case of H = 0.3 (where the crossing point 
moves rightward). 

Clearly in this model (and probably in any spin glass 
model in a field), there are finite size effects with oppo- 
site signs in competition, whose relative weights depend 
on the value of the external field. This makes the extrap- 
olation to the thermodynamic limit extremely challeng- 
ing, especially in cases like the one in Fig. |4] (middle) 
where any reasonable extrapolation would predict a crit- 
ical temperature (if any) much lower than the true one. 
This phenomenon must be taken into account when one 
wants to exclude any phase transition because of lack of 
a crossi ng p oint in the scaled data. See, for example, 
Refs. I25ll26l for a recent discussion about this issue for 
mean-field and non-mean-field spin glasses in an exter- 
nal field. 

Despite the large finite size effects, we have tried to 
scale the data according to Eg . (|20|) . The results are ob- 
viously very poor and have been reported in Fig.[5]for the 
sake of completeness. Even optimizing over the choice of 
the 7 parameter (such as to superimpose the data from 
largest sizes at the critical point) we get no data collapse 
at all (sec main panels in Fig. [5]). The tentative scaling 
with 7 = 1/2 is shown in the insets of Fig. [5] and it is 
even worst. 
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FIG. 4: Rescaled spin glass susceptibility (assuming 7 = 1/2) 
versus temperature for fields H — 0.1 (top), H = 0.2 (middle), 
and H = 0.3 (bottom). Errors are smaller than the symbol 
size. 



We should note en passant that a much better scal- 
ing of the numerical data can be obtained by using two 
different exponents for the rescalings of the x- and y- 
axes. This scaling would imply a divergence of the SG 
susceptibility in a field as xsg ^ \t\~ a with a ~ 1.6 -r 1.7. 
However, given the strong analytical arguments in favor 
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FIG. 5: "Tentative" scaling function for the spin glass sus- 
ceptibility for fields H — 0.1 (top), H — 0.2 (middle), and 
H = 0.3 (bottom). In the main panels, the parameter 7 has 
been optimized such as to superimpose the data from the 
largest sizes at the critical point. In the insets, the value 
7 = 1/2 is used. Errors are smaller than the symbol size. 
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H 





0.1 


0.2 


0.3 


7 (Ai) 


0.598±0.012 


0.545±0.010 


0.943±0.014 


0.789±0.015 


7 (A?) 


0.599±0.013 


0.589±0.021 


1.139±0.001 


0.730±0.033 



TABLE II: Values of 7 estimated from the arithmetic averages 
over 2000 samples of Ai (Ai) and those of \\ (Af). Errors 
represent the standard errors of the estimates. For H = 0, 
the two estimates coincide with each other up to the second 
digit, justifying a scaling of the form Ai = iV _1 ^' 1+7 ^i, where 
£1 is a certain random variable that is independent of N. On 
the other hand, they differ significantly for H = 0.1, 0.2, and 
0.3, which implies that a scaling of that form does not hold 
for H > 0. 

of a = 1, we have to conclude that such an alternative 
scaling is only due to finite size effects. 

In order to examine the reason for the poor consis- 
tency with the scaling law with 7 = 1/2, i.e. xsg = 
A rl / 3 F(^iV 1 / 3 ), we numerically evaluated the eigenval- 
ues of the Hessian matrix A = for systems of 
N = 2 5 ,2 6 ,2 7 , and 2 8 . The reduction of the system 
size is simply due to limited computational resources; the 
eigenvalue analysis costs much more than the evaluation 
of Xsg- 

Figures [H] (a)-(d) show the cumulative distributions 
Jq dtp(t) at the critical temperature T c , which were 
numerically computed from data of 2000 sample sys- 
tems. These distributions should be proportional to 
(A — A m ; n ) 1+7 — ► A 1+7 in the limit of TV —5- 00 since 
Amin vanishes at criticality. However, finite size correc- 
tions modify the behavior of the numerical data in the 
vicinity of A = 0. The envelopes of the cumulative dis- 
tributions for various TV exhibit reasonable consistency 
with the scaling form of 7 = 1/2, which is represented by 
straight lines in Figs. [HI for all cases of H = 0, 0.1, 0.2, 
and 0.3. 

Another possible source of inconsistency concerning 
the finite size scaling relation of xsg is the finite size 
correction to the lower band edge A m ; n . The argu- 
ment presented in Sec. IIII Bl relies on the assumption 
that the smallest eigenvalue Ai behaves as Ai ~ A min + 
~ at T c , where £1 is a random 

variable independent of N. This means that the distri- 
butions of Ai for different A" collapse to a single curve 
for fixed H after rescaling of Ai -> A^/^+^Ai using an 
appropriate value of 7. Figure [7] (a) shows that this holds 
to good accuracy for H = with 7 = 0.598, which is rea- 
sonably close to the theoretical prediction 1/2. However, 
Figs. [7] (b)-(d) and Table |TT] indicate that such a clear 
scaling relation does not hold for H > 0. 

This is presumably because, in the presence of an ex- 
ternal field, the lower band edge A lllm has strong finite 
size corrections and a non-negligible dependence on N. 
The critical temperature T c is defined by the condition 
Amin = in the limit of N —> 00, but for finite N in 
the presence of an external field, it may not vanish due 
to a positive bias as A min = B 1 N~ a even at T c . Such 
a bias could come out if statistical correlations among 
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FIG. 6: Cumulative eigenvalue distribution dtp(t) at the 
critical temperature T c for (a): H — 0, (b): H — 0.1, (c): 
H = 0.2, and (d): H = 0.3. From right to left, the curves 
correspond to N = 2 5 , 2 6 , 2 7 , and 2 8 . The straight lines stand 
for a scaling relation of dtp(t) oc A 1+7 with 7 = 1/2. 



FIG. 7: Scaling plots for the cumulative distributions of the 
smallest eigenvalue Ai at the critical temperatures T c for (a): 
H = 0, (b): H = 0.1, (c): H = 0.2, and (d): H = 0.3. In 
each plot, four curves correspond to N — 2 5 , 2 6 , 2 7 , and 2 s . 
The plots are obtained by rescaling the horizontal axis of the 
original ones (insets: N = 2 5 , 2 6 , 2 7 , and 2 s from right to 
left) as Ai — > iV 1/ '' 1+7 'Ai, where 1/(1 +7) is determined on 
the basis of the arithmetic averages of Ai over 2000 samples. 
The estimates of 7 are 0.598, 0.545, 0.943, and 0.789 for 
H = 0.0, 0.1, 0.2, and 0.3, respectively. Except for (a), the 
data do not collapse to a single curve with good accuracy. 
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the eigenvalues Ai, Aa, ■ ■ ■ , Xn are not negligible. Con- 
sequently, the expression for the smallest eigenvalue Ai 
at T c should be modified at least as 

Ai = BiN-o +£i7V~ 1/(1+7) , 

where random variable £1 generally obeys a certain non- 
trivial distribution whose mean is not guaranteed to 
vanish^. Unless a and 1/(1 + 7) arc very different (or 
very similar), both scaling terms will be needed for appro- 
priately handling data with N of several hundreds, which 
are the practical upper limits of the system size that we 
can deal with by standard computational resources to 
date. However, estimating the two exponents simultane- 
ously from data of only few values of N, is far from trivial 
in the presence of statistical fluctuations. 

The presence of strong finite size corrections for A m i n 
and the practical difficulty of identifying the scaling re- 
lation of the corrections from numerical data mean that 
accurate numerical evaluation of T c is very difficult in the 
presence of external fields. Although we herein examined 
random sparse networks, a similar issue should also affect 
other systems. This may be a major reason why the AT 
line has not been clearly observed in SG models of finite 
dimensions. 



V. SUMMARY 

In summary, we have explored a finite size scaling re- 
lation of the de Almeida-Thouless (AT) instability crit- 
icality in random sparse networks analytically and nu- 
merically. The spin glass susceptibility xsg is a proper 
measure for signaling the AT criticality. On the basis of 
the similarity between the random sparse networks and 
the Bcthe trees, we have derived a scheme for evaluating 
Xsg of infinitely large systems utilizing the belief prop- 
agation algorithm, which makes it possible to evaluate 
the critical temperature T c with a numerically feasible 
procedure. The singularity at T c , however, cannot be di- 
rectly observed in finite systems due to finite size effects. 
Therefore, we examined how the finite size scaling rela- 
tion is determined by the lower band edge behavior of 
the Hessian matrix. 

The validity of the theoretical predictions was ex- 
amined in extensive numerical experiments. For suffi- 
ciently high temperatures, the numerically computed val- 
ues of xsg were reasonably consistent with the theoreti- 
cal predictions regardless of whether an external field was 
present. This result supports a widely believed equiva- 
lence between random sparse networks and Bcthe trees, 
implying that T c is the same in both systems. Accord- 
ingly, we investigated the consistency between the numer- 
ical data and the theoretically obtained finite size scaling 
relation XSG = N w F (N^ \t\) . In the absence of an exter- 
nal field, the numerical data are in good agreement with 
the scaling relation of u = 1/3, as was believed in earlier 
studies. On the other hand, the consistency of the data 
with respect to the scaling relation becomes very poor in 
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FIG. 8: Rescaled spin glass susceptibility (assuming 7 = 1/2 
as in Fig. [4j versus temperature. The critical line T C {H) is 
crossed perpendicularly. The finite size effects are smaller 
than in Fig. [?] (middle), but have the same qualitative behav- 
ior. 

the presence of a field; the crossing points of -/V _1//3 xsg 
for N < 2 10 fluctuate around the theoretical values of T c 
non-monotonically with the strength of the field H, and 
the data do not fit the scaling relation even if w is tuned. 
Upon examining the eigenvalues of the Hessian matrix on 
the basis of numerical simulations for N < 2 s , we found 
that the lower band edge of the eigenvalue distribution 
is very sensitive to H at least for N of several hundreds. 
This might be a major reason for the inconsistency with 
the theoretical prediction of finite size scaling. 

Reference suggested studying the AT instability 
along a path in the (T, H) plane that perpendicularly 
crosses the critical line T C (H), which successfully led to 
accurate estimates of the criticality based on data of 
N < 2 9 in the case of C = 6: in this way, the finite 
size effect should be reduced. We followed such a sugges- 
tion and analyzed our data along the path represented 
by the dashed line in the inset of Fig. [HJ which crosses 
the critical line (full line in the inset of Fig. [8j perpen- 
dicularly at H = 0.2. The resulting susceptibility, scaled 
by the factor iV -1 / 3 as in Fig. 21 is shown in the main 
panel of Fig. [8J and should be compared with the data 
reported in Fig. 2] (middle). It is interesting to note that 
finite size effects are indeed much reduced (roughly by a 
factor 4), but the qualitative behavior of the data is ex- 
actly the same as in the case with H fixed. In particular, 
the crossing temperature moves away from the critical 
temperature, thus giving too small a T c . Certainly for 
larger sizes, the crossing temperature will come back to 
T c but we have no numerical evidence of that for sizes 
up to N — 2 10 . The same behavior occurs if the data is 
plotted against the field intensity, as in Fig. 4 of Ref. [^. 

In light of the present results on the existence of very 
strong finite size effects, we believe that the outcome of 
numerical simulations of spin glasses in a field should be 
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taken with a lot of care. In particular, the observation 
that rescaled susceptibilities do not cross in a wide tem- 
perature ranged or have a crossing point moving to low 
temperatures^ should not be taken as a definite indica- 
tion for the lack of a spin glass phase. In the present 
work we have shown that even in mean field models very 
strong and possibly non-monotonic corrections to finite- 
size scaling exist. In finite dimensional models these cor- 
rections are likely to become stronger and extrapolation 
from relatively small system sizes is risky. Maybe the de- 
velopment of new data analysis methods that may help 



to reduce finite size effects (like the one in Ref. |26| ) would 
be very welcome. 
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